/*  Quasi 1 - D Euler Equation Solver for Laval Tunnel by FVM
 */
#include <iostream>
#include <vector>
#include <cmath>
#include <Eigen/Dense>
#include <eigen3/unsupported/Eigen/CXX11/Tensor>
using namespace std;
using namespace Eigen;

class laval
{
public:
    enum
    {
        RK3,
        RK4
    };

    int ncells_;                  // num of cells
    double dx_, xl_, xr_, a_, dt; // domain and dx
    double cfl;
    constexpr static double gamma = 1.4;
    constexpr static double k2 = 0.75;
    constexpr static double k4 = 1.0 / 64.0;
    // constexpr static double cfl = 0.85;
    // constexpr static double dt = 0.001;
    //  ----------------- bc & far & lav.ncells_ + 1 -------------
    constexpr static double r0 = 1.0;
    constexpr static double p0 = 1.0, pe = 0.8785;

    // ---------------------- flux ---------------------
    double (*U_)[3]; // U
    double (*F_)[3]; // F
    double (*H_)[3]; // H
    double (*R_)[3]; // res
    double (*D_)[3]; // vis

    // ---------------------- var ---------------------
    vector<double> xc_;   // cordination of cell center
    vector<double> u_;    // speed
    vector<double> r_;    // density
    vector<double> p_;    // pressure
    vector<double> e_;    // energy
    vector<double> m_;    // ma
    vector<double> q_;    // sound speed
    vector<double> mass_; // mass flux

    // ---------------------- vis ---------------------
    vector<double> alpha_;
    vector<double> e2_, e4_;
    vector<double> theta_;

public:
    laval(double xl, double xr, int ncells) : xl_(xl), xr_(xr), dx_((xr - xl) / ncells), ncells_(ncells),
                                              U_(new double[ncells][3]),
                                              F_(new double[ncells][3]),
                                              H_(new double[ncells][3]),
                                              R_(new double[ncells][3]),
                                              D_(new double[ncells][3]),
                                              u_(ncells + 2), r_(ncells + 2), p_(ncells + 2), e_(ncells + 2),
                                              e2_(ncells + 1), e4_(ncells + 1), theta_(ncells), mass_(ncells),
                                              xc_(ncells + 2), m_(ncells + 2), q_(ncells + 2), alpha_(ncells + 1)

    {
        // ---------------------- func ---------------------
    }

public:
    void cal_e(void) // calculate energy
    {
        for (int i = 1; i < ncells_ + 1; i++)
            e_[i] = 1.0 / (gamma - 1) * (p_[i] / r_[i]) + 0.5 * u_[i] * u_[i];
        e_[0] = 1.0 / (gamma - 1) * (p_[0] / r_[0]) + 0.5 * u_[0] * u_[0];
        e_[ncells_ + 1] = 1.0 / (gamma - 1) * (p_[ncells_ + 1] / r_[ncells_ + 1]) + 0.5 * u_[ncells_ + 1] * u_[ncells_ + 1];
        return;
    }

    /**
     * @brief output in .plt file
     */
    void tecout(void)
    {
        FILE *fp = fopen("Laval.plt", "w");

        fprintf(fp, "variables = \"x\", \"u\"\n");
        fprintf(fp, "Zone T = \"u\"\n");
        for (int i = 1; i < ncells_ + 2; i++)
            fprintf(fp, "%lf  %lf\n", xc_[i], u_[i]);

        fprintf(fp, "variables = \"x\", \"p\"\n");
        fprintf(fp, "Zone T = \"p\"\n");
        for (int i = 1; i < ncells_ + 2; i++)
            fprintf(fp, "%lf  %lf\n", xc_[i], p_[i]);

        fprintf(fp, "variables = \"x\", \"Ma\"\n");
        fprintf(fp, "Zone T = \"Ma\"\n");
        for (int i = 1; i < ncells_ + 2; i++)
            fprintf(fp, "%lf  %lf\n", xc_[i], m_[i]);

        fprintf(fp, "variables = \"x\", \"Mass\"\n");
        fprintf(fp, "Zone T = \"Mass\"\n");
        for (int i = 1; i < ncells_; i++)
            fprintf(fp, "%lf  %lf\n", xc_[i], mass_[i]);

        fp = NULL;
        return;
    }

    /**
     * @brief calculate flux
     * @param rho
     * @param u
     * @param p
     * @return flux
     */
    void flux(double (*U)[3])
    {
        //< var for interface
        VectorXd uu(ncells_ + 1);
        VectorXd rr(ncells_ + 1);
        VectorXd pp(ncells_ + 1);
        VectorXd ee(ncells_ + 1);

        for (int i = 1; i < ncells_; i++)
        {
            uu[i] = 0.5 * (u_[i + 1] + u_[i]);
            rr[i] = 0.5 * (r_[i + 1] + r_[i]);
            pp[i] = 0.5 * (p_[i + 1] + p_[i]);
            ee[i] = 1 / (gamma - 1) * (p_[i] / r_[i]) + 0.5 * u_[i] * u_[i];
        }
        uu[0] = (u_[0] + u_[1]) / 2.0; //< (u_inf + u_1)*0.5
        rr[0] = (r_[0] + r_[1]) / 2.0;
        pp[0] = (p_[0] + p_[1]) / 2.0;
        // uu[0] = u_[0];
        // rr[0] = r_[0];
        // pp[0] = p_[0];
        ee[0] = 1.0 / (gamma - 1) * (p_[0] / r_[0]) + 0.5 * u_[0] * u_[0];

        uu[ncells_] = u_[ncells_ + 1];
        rr[ncells_] = r_[ncells_ + 1];
        pp[ncells_] = p_[ncells_ + 1];
        ee[ncells_] = 1.0 / (gamma - 1) * (p_[ncells_ + 1] / r_[ncells_ + 1]) + 0.5 * u_[ncells_ + 1] * u_[ncells_ + 1];

        for (int i = 1; i < ncells_ - 1; i++)
        {
            U[i][0] = r_[i + 1];
            U[i][1] = r_[i + 1] * u_[i + 1];
            U[i][2] = r_[i + 1] * e_[i + 1];

            // std::cout << uu[i] << endl;
            F_[i][0] = rr[i + 1] * uu[i + 1] - rr[i] * uu[i]; // calculate flux F
            F_[i][1] = rr[i + 1] * uu[i + 1] * uu[i + 1] + pp[i + 1] -
                       rr[i] * uu[i] * uu[i] - pp[i];
            F_[i][2] = (rr[i + 1] * ee[i + 1] + pp[i + 1]) * uu[i + 1] -
                       (rr[i] * ee[i] + pp[i]) * uu[i];

            H_[i][0] = r_[i + 1] * u_[i + 1] * 2.0 * xc_[i + 1] / (0.5 + xc_[i + 1] * xc_[i + 1]) * dx_;
            H_[i][1] = r_[i + 1] * u_[i + 1] * u_[i + 1] * 2.0 * xc_[i + 1] / (0.5 + xc_[i + 1] * xc_[i + 1]) * dx_;
            H_[i][2] = (r_[i + 1] * e_[i + 1] + p_[i + 1]) * u_[i + 1] * 2.0 * xc_[i + 1] / (0.5 + xc_[i + 1] * xc_[i + 1]) * dx_;
        }

        //< special work for first & last cells
        U[0][0] = r_[1];
        U[0][1] = r_[1] * u_[1];
        U[0][2] = r_[1] * e_[1];

        F_[0][0] = rr[1] * uu[1] - rr[0] * uu[0];
        F_[0][1] = rr[1] * uu[1] * uu[1] + pp[1] -
                   rr[0] * uu[0] * uu[0] - pp[0];
        F_[0][2] = (rr[1] * ee[1] + pp[1]) * uu[1] -
                   (rr[0] * ee[0] + pp[0]) * uu[0];

        H_[0][0] = r_[1] * u_[1] * 2.0 * xc_[1] / (0.5 + xc_[1] * xc_[1]) * dx_;
        H_[0][1] = r_[1] * u_[1] * u_[1] * 2.0 * xc_[1] / (0.5 + xc_[1] * xc_[1]) * dx_;
        H_[0][2] = (r_[1] * e_[1] + p_[1]) * u_[1] * 2.0 * xc_[1] / (0.5 + xc_[1] * xc_[1]) * dx_;

        U[ncells_ - 1][0] = r_[ncells_];
        U[ncells_ - 1][1] = r_[ncells_] * u_[ncells_];
        U[ncells_ - 1][2] = r_[ncells_] * e_[ncells_];

        F_[ncells_ - 1][0] = rr[ncells_] * uu[ncells_] - rr[ncells_ - 1] * uu[ncells_ - 1];
        F_[ncells_ - 1][1] = rr[ncells_] * uu[ncells_] * uu[ncells_] + pp[ncells_] -
                             rr[ncells_ - 1] * uu[ncells_ - 1] * uu[ncells_ - 1] - pp[ncells_ - 1];
        F_[ncells_][2] = (rr[ncells_] * ee[ncells_] + pp[ncells_]) * uu[ncells_] -
                         (rr[ncells_ - 1] * ee[ncells_ - 1] + pp[ncells_ - 1]) * uu[ncells_ - 1];

        H_[ncells_ - 1][0] = r_[ncells_] * u_[ncells_] * 2.0 * xc_[ncells_] / (0.5 + xc_[ncells_] * xc_[ncells_]) * dx_;
        H_[ncells_ - 1][1] = r_[ncells_] * u_[ncells_] * u_[ncells_] * 2.0 * xc_[ncells_] / (0.5 + xc_[ncells_] * xc_[ncells_]) * dx_;
        H_[ncells_ - 1][2] = (r_[ncells_] * e_[ncells_] + p_[ncells_]) * u_[ncells_] * 2.0 * xc_[ncells_] / (0.5 + xc_[ncells_] * xc_[ncells_]) * dx_;

        return;
    }

    /**
     * @brief calculate res for time advance
     * @param Flux
     * @param Vis
     * @return Res
     */
    void res(void)
    {
        MatrixXd d(ncells_ + 1, 3); //< vis for interface
        MatrixXd U_g(3, 3);         //< ghost cells
        double theta_b1[3];         //< theta[-1], theta[i+1], theta[i+2]

        //< ---------------------------------------------------special work for vis in first & last two interface
        alpha_[0] = 0.5 * (dx_ / dt + dx_ / dt); // ------------------------------
        theta_[0] = fabs((p_[2] - 2.0 * p_[1] + p_[0]) / (p_[2] + 2.0 * p_[1] + p_[0]));
        theta_b1[0] = fabs((p_[1] - 2.0 * p_[0] + p_[0]) / (p_[1] + 2.0 * p_[0] + p_[0]));
        U_g(0, 0) = r_[0];
        U_g(0, 1) = r_[0] * u_[0];
        U_g(0, 2) = r_[0] * 1.0 / (gamma - 1) * (p_[0] / r_[0]) + 0.5 * u_[0] * u_[0];
        e2_[0] = k2 * max(max(theta_b1[0], theta_[0]), theta_[1]);
        e4_[0] = max(0.0, k4 - e2_[0]); //--------------------------------------
        for (int j = 0; j < 3; j++)     //< first cell
            d(0, j) = alpha_[0] * (e2_[0] * (U_[0][j] - U_g(0, j)) - e4_[0] * (U_[1][j] - 3.0 * U_[0][j] + 3.0 * U_g(0, j) - U_g(0, j)));

        alpha_[ncells_] = 0.5 * (dx_ / dt + dx_ / dt); // ------------------------------
        theta_b1[1] = fabs((p_[ncells_ + 1] - 2.0 * p_[ncells_ + 1] + p_[ncells_]) / (p_[ncells_ + 1] + 2.0 * p_[ncells_ + 1] + p_[ncells_]));
        theta_b1[2] = 0.0;
        theta_[ncells_ - 1] = fabs((p_[ncells_ + 1] - 2.0 * p_[ncells_] + p_[ncells_ - 1]) / (p_[ncells_ + 1] + 2.0 * p_[ncells_] + p_[ncells_ - 1]));
        U_g(1, 0) = r_[ncells_ + 1];
        U_g(1, 1) = r_[ncells_ + 1] * u_[ncells_ + 1];
        U_g(1, 2) = r_[ncells_ + 1] * 1.0 / (gamma - 1) * (p_[ncells_ + 1] / r_[ncells_ + 1]) + 0.5 * u_[ncells_ + 1] * u_[ncells_ + 1];
        U_g(2, 0) = r_[ncells_ + 1];
        U_g(2, 1) = r_[ncells_ + 1] * u_[ncells_ + 1];
        U_g(2, 2) = r_[ncells_ + 1] * 1.0 / (gamma - 1) * (p_[ncells_ + 1] / r_[ncells_ + 1]) + 0.5 * u_[ncells_ + 1] * u_[ncells_ + 1];
        e2_[ncells_] = k2 * max(max(max(theta_[ncells_ - 2], theta_[ncells_ - 1]), theta_b1[1]), theta_b1[2]);
        e4_[ncells_] = max(0.0, k4 - e2_[ncells_]);

        alpha_[ncells_ - 1] = 0.5 * (dx_ / dt + dx_ / dt); // ------------------------------------------
        theta_[ncells_ - 2] = fabs((p_[ncells_] - 2.0 * p_[ncells_ - 1] + p_[ncells_ - 2]) / (p_[ncells_] + 2.0 * p_[ncells_ - 1] + p_[ncells_ - 2]));
        e2_[ncells_ - 1] = k2 * max(max(max(theta_[ncells_ - 3], theta_[ncells_ - 2]), theta_[ncells_ - 1]), theta_b1[1]); // last tow cell
        e4_[ncells_ - 1] = max(0.0, k4 - e2_[ncells_ - 1]);

        for (int j = 0; j < 3; j++) // last two cell
        {
            d(ncells_ - 1, j) = alpha_[ncells_ - 1] *
                                (e2_[ncells_ - 1] * (U_[ncells_ - 1][j] - U_[ncells_ - 2][j]) - e4_[ncells_ - 1] * (U_g(1, j) - 3.0 * U_[ncells_ - 1][j] + 3.0 * U_[ncells_ - 2][j] - U_[ncells_ - 3][j]));
        }

        for (int j = 0; j < 3; j++) // last cell
        {
            d(ncells_, j) = alpha_[ncells_] *
                            (e2_[ncells_] * (U_g(1, j) - U_[ncells_ - 1][j]) - e4_[ncells_] * (U_g(2, j) - 3.0 * U_g(1, j) + 3.0 * U_[ncells_ - 1][j] - U_[ncells_ - 2][j]));
        }

        //< ---------------------------------------- inner cells
        for (int i = 2; i < ncells_ - 1; i++)
        {
            theta_[i] = fabs((p_[i + 2] - 2.0 * p_[i + 1] + p_[i]) / (p_[i + 2] + 2.0 * p_[i + 1] + p_[i]));
            alpha_[i] = 0.5 * (dx_ / dt + dx_ / dt);
            e2_[i] = k2 * max(max(max(theta_[i - 2], theta_[i - 1]), theta_[i]), theta_[i + 1]);
            e4_[i] = max(0.0, k4 - e2_[i]);
            for (int j = 0; j < 3; j++)
                d(i, j) = alpha_[i] * (e2_[i] * (U_[i][j] - U_[i - 1][j]) - e4_[i] * (U_[i + 1][j] - 3.0 * U_[i][j] + 3.0 * U_[i - 1][j] - U_[i - 2][j]));
        }

        //-----------------------------------first two cell
        alpha_[1] = 0.5 * (dx_ / dt + dx_ / dt); // ------------------------------
        e2_[1] = k2 * max(max(max(theta_b1[0], theta_[0]), theta_[1]), theta_[2]);
        theta_[1] = fabs((p_[3] - 2.0 * p_[2] + p_[1]) / (p_[3] + 2.0 * p_[2] + p_[1]));
        for (int j = 0; j < 3; j++) //< first 2 cell
            d(1, j) = alpha_[1] * (e2_[1] * (U_[1][j] - U_[0][j]) - e4_[1] * (U_[2][j] - 3.0 * U_[1][j] + 3.0 * U_[0][j] - U_g(0, j)));

        //__________________________________________________________________________________________

        //< calculate res
        for (int i = 0; i < ncells_; i++)
        {
            for (int j = 0; j < 3; j++)
            {
                D_[i][j] = d(i + 1, j) - d(i, j);
                R_[i][j] = (F_[i][j] + H_[i][j] - D_[i][j]) / dx_;
            }
        }

        return;
    }

#if 0
    /**
 * @brief U to rho, u, e; update in each step 
 * @param U
 * @return var
 */
    void u2var_sub(double (*U)[3])
    {
        for (int i = 1; i < ncells_ + 1; i++)
        {
            r_[i] = U[i - 1][0];
            u_[i] = U[i - 1][1] / r_[i];
            e_[i] = U[i - 1][2] / r_[i];
            p_[i] = r_[i] * (gamma - 1) * (e_[i] - 0.5 * u_[i] * u_[i]);
        }
        u_[0] = (3 * u_[1] - u_[2]) * 0.5; // ------------------- update bc cells
        u_[ncells_ + 1] = (3 * u_[ncells_] - u_[ncells_ - 1]) * 0.5;
        r_[0] = pow(0.4 / 1.4 * (1.4 / 0.4 - 0.5 * u_[0] * u_[0]), 2.5);
        p_[0] = pow(r_[0], 1.4);
        p_[ncells_ + 1] = 0.8785;
        r_[ncells_ + 1] = pow(r_[ncells_], 1.4) / p_[ncells_] * p_[ncells_ + 1];
        r_[ncells_ + 1] = pow(r_[ncells_ + 1], 2.5);
        q_[ncells_ + 1] = sqrt(gamma * p_[ncells_ + 1] / r_[ncells_ + 1]);
        m_[ncells_ + 1] = u_[ncells_ + 1] / q_[ncells_ + 1];
    }
#endif

    /**
     * @brief 用于在时间推进的每一个子步计算原始变量(入口超声速，出口亚声速)
     * @param 守恒变量U
     * @return 密度、压力、速度
     */
    void u2var_super(double (*U)[3])
    {
        for (int i = 1; i < ncells_ + 1; i++)
        {
            r_[i] = U[i - 1][0];
            u_[i] = U[i - 1][1] / r_[i];
            e_[i] = U[i - 1][2] / r_[i];
            p_[i] = r_[i] * (gamma - 1) * (e_[i] - 0.5 * u_[i] * u_[i]);
            q_[i] = sqrt(gamma * p_[i] / r_[i]);
            m_[i] = u_[i] / q_[i];
        }
        m_[0] = 1.16;
        p_[0] = 1.0 / pow((1.0 + 0.2 * m_[1] * m_[1]), 3.5);
        r_[0] = 1.0 / pow((1.0 / p_[1]), 1.0 / 1.4);
        u_[0] = m_[0] * sqrt(1.4 * p_[0] / r_[0]);

        // p_[0] = 0.5 * (p_[1] + 1.0);
        // r_[0] = 0.5 * (r_[1] + 1.0);
        // u_[0] = 0.5 * (u_[1] + u_[0]);

        p_[ncells_ + 1] = 0.8785;
        u_[ncells_ + 1] = (3 * u_[ncells_] - u_[ncells_ - 1]) * 0.5;
        r_[ncells_ + 1] = pow(r_[ncells_], 1.4) * p_[ncells_ + 1] / p_[ncells_];
        r_[ncells_ + 1] = pow(r_[ncells_ + 1], 1 / 1.4);
        q_[ncells_ + 1] = sqrt(gamma * p_[ncells_ + 1] / r_[ncells_ + 1]);
        m_[ncells_ + 1] = u_[ncells_ + 1] / q_[ncells_ + 1];

        return;
    }

    /**
     * @brief time advance by rk4
     */
    void rk4(void)
    {
        double(*U_temp)[3];
        Tensor<double, 3> UU(ncells_, 3, 5);
        vector<double>::iterator loc_u;
        U_temp = new double[ncells_][3];

        // cal_e();
        flux(U_);
        res();

        for (int i = 0; i < ncells_; i++)
        {
            for (int j = 0; j < 3; j++)
            {
                UU(i, j, 0) = U_[i][j];
                UU(i, j, 1) = UU(i, j, 0) - 0.25 * dt * R_[i][j];
                U_temp[i][j] = UU(i, j, 1); //< use temp here to store
            }
        }
        u2var_super(U_temp);
        flux(U_temp);
        res(); //< update res

        for (int i = 0; i < ncells_; i++)
        {
            for (int j = 0; j < 3; j++)
            {
                UU(i, j, 2) = U_[i][j] - 1.0 / 3.0 * dt * R_[i][j];
                U_temp[i][j] = UU(i, j, 2);
            }
        }
        u2var_super(U_temp);
        flux(U_temp);
        res(); //< update res

        for (int i = 0; i < ncells_; i++)
        {
            for (int j = 0; j < 3; j++)
            {
                UU(i, j, 3) = U_[i][j] - 0.5 * dt * R_[i][j];
                U_temp[i][j] = UU(i, j, 3);
            }
        }
        u2var_super(U_temp);
        flux(U_temp);
        res(); //< update res

        for (int i = 0; i < ncells_; i++)
        {
            for (int j = 0; j < 3; j++)
            {
                UU(i, j, 4) = U_[i][j] - dt * R_[i][j];
                U_[i][j] = UU(i, j, 4); //< update var
            }
            mass_[i] = r_[i + 1] * u_[i + 1] * (0.5 + xc_[i + 1] * xc_[i + 1]);
        }
        delete[] U_temp;
        u2var_super(U_);
        flux(U_);
        res(); //< update res

        loc_u = max_element(u_.begin(), u_.end());

        dt = cfl * dx_ / (1.5 * *loc_u);

        return;
    }

    /**
     * @brief time advance RK3
     */
    void rk3(void)
    {
        double(*U_temp)[3];
        Tensor<double, 3> UU(ncells_, 3, 4);
        vector<double>::iterator loc_u;
        U_temp = new double[ncells_][3];

        flux(U_);
        res();

        for (int i = 0; i < ncells_; i++)
        {
            for (int j = 0; j < 3; j++)
            {
                UU(i, j, 0) = U_[i][j];
                UU(i, j, 1) = UU(i, j, 0) - dt * R_[i][j];
                U_temp[i][j] = UU(i, j, 1);
            }
        }

        u2var_super(U_temp);
        flux(U_temp);
        res(); //< update res

        for (int i = 0; i < ncells_; i++)
        {
            for (int j = 0; j < 3; j++)
            {
                UU(i, j, 2) = 0.75 * UU(i, j, 0) + 0.25 * UU(i, j, 1) - 0.25 * dt * R_[i][j];
                U_temp[i][j] = UU(i, j, 2);
            }
        }
        u2var_super(U_temp);
        flux(U_temp);
        res(); //< update res

        for (int i = 0; i < ncells_; i++)
        {
            for (int j = 0; j < 3; j++)
            {
                UU(i, j, 3) = 1.0 / 3.0 * UU(i, j, 0) + 2.0 / 3.0 * UU(i, j, 2) - 2.0 / 3.0 * dt * R_[i][j];
                U_[i][j] = UU(i, j, 3);
            }
            mass_[i] = r_[i + 1] * u_[i + 1] * (0.5 + xc_[i + 1] * xc_[i + 1]);
        }
        delete[] U_temp;
        u2var_super(U_);
        flux(U_);
        res(); //< update res

        loc_u = max_element(u_.begin(), u_.end());

        dt = cfl * dx_ / (*loc_u);

        return;
    }

    void time_advance(int method)
    {
        if (method == RK3)
            rk3();
        else if (method == RK4)
            rk4();
    }
};

#if 0
/** 
 * @brief 入口出口均为亚声速的边界条件
 */
void sub_init(laval &lav, vector<double> &x_, vector<double> &r_, vector<double> &p_, vector<double> &e_,
              vector<double> &u_, vector<double> &m_, vector<double> &q_)
{
    vector<double> t(lav.ncells_ + 2);
    for (int i = 1; i < lav.ncells_ + 1; i++)
    {
        x_[i] = lav.xl_ + (i - 1) * lav.dx_ + 0.5 * lav.dx_;

        p_[i] = 1.0 + (i - 1) * (1.0 - 0.8785) / (lav.ncells_);
        r_[i] = pow(p_[i], 1.0 / 1.4);
        t[i] = 1.0 - 0.2314 * x_[i];
        u_[i] = sqrt(fabs(2.0 * (3.5 - 3.5 * (p_[i] / r_[i]))));

        e_[i] = 1.0 / (lav.gamma - 1) * (p_[i] / r_[i]) + 0.5 * u_[i] * u_[i];
        q_[i] = sqrt(lav.gamma * p_[i] / r_[i]);
        m_[i] = u_[i] / q_[i];
    }

    x_[0] = -1.0 / 3.0;
    x_[lav.ncells_ + 1] = 1.0;
    lav.dt = 0.001;
    lav.cfl = 0.85;

    u_[0] = (3 * u_[1] - u_[2]) * 0.5; // ------------------- update in iteration !
    u_[lav.ncells_ + 1] = (3 * u_[lav.ncells_] - u_[lav.ncells_ - 1]) * 0.5;
    r_[0] = pow((1.0 - 0.4 / 1.4 * 0.5 * u_[0] * u_[0]), 2.5);
    p_[0] = pow(r_[0], 1.4);
    p_[lav.ncells_ + 1] = 0.8785; //(3 * p_[lav.ncells_] - p_[lav.ncells_ - 1]) * 0.5;
    r_[lav.ncells_ + 1] = pow(r_[lav.ncells_], 1.4) * p_[lav.ncells_ + 1] / p_[lav.ncells_];
    r_[lav.ncells_ + 1] = pow(r_[lav.ncells_ + 1], 1 / 1.4);

    return;
}
#endif

/**
 * @brief 入口超声速，出口亚声速的边界条件设置
 */
void super_init(laval &lav, vector<double> &x_, vector<double> &r_, vector<double> &p_, vector<double> &e_,
                vector<double> &u_, vector<double> &m_, vector<double> &q_)
{
    vector<double> t(lav.ncells_ + 2);

    for (int i = 1; i < lav.ncells_ + 2; i++)
    {
        x_[i] = lav.xl_ + (i - 1) * lav.dx_ + 0.5 * lav.dx_;

        p_[i] = 1.0 + (i - 1) * (1.0 - 0.8785) / (lav.ncells_);
        r_[i] = pow(p_[i], 1.0 / 1.4);
        t[i] = 1.0 - 0.2314 * x_[i];
        u_[i] = sqrt(fabs(2.0 * (3.5 - 3.5 * (p_[i] / r_[i]))));

        e_[i] = 1.0 / (lav.gamma - 1) * (p_[i] / r_[i]) + 0.5 * u_[i] * u_[i];
        q_[i] = sqrt(lav.gamma * p_[i] / r_[i]);
        m_[i] = u_[i] / q_[i];
    }

    x_[0] = 0.1;
    x_[lav.ncells_ + 1] = 1.00;
    lav.dt = 0.001;
    lav.cfl = 0.6;

    m_[0] = 1.16;
    p_[0] = 1.0 / pow((1.0 + 0.2 * m_[1] * m_[1]), 3.5);
    r_[0] = 1.0 / pow((1.0 / p_[1]), 1.0 / 1.4);
    u_[0] = m_[0] * sqrt(1.4 * p_[0] / r_[0]);
    q_[0] = sqrt(lav.gamma * p_[0] / r_[0]);

    p_[lav.ncells_ + 1] = 0.8785;
    u_[lav.ncells_ + 1] = (3 * u_[lav.ncells_] - u_[lav.ncells_ - 1]) * 0.5;
    r_[lav.ncells_ + 1] = pow(r_[lav.ncells_], 1.4) * p_[lav.ncells_ + 1] / p_[lav.ncells_];
    r_[lav.ncells_ + 1] = pow(r_[lav.ncells_ + 1], 1 / 1.4);
    q_[lav.ncells_ + 1] = sqrt(lav.gamma * p_[lav.ncells_ + 1] / r_[lav.ncells_ + 1]);

    return;
}